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f*****"**"*^®®^»^**^  Prof .  R.  Krishnamurti  and  Dr  H.  Yang  have  completed  their 
investigation  of  finite  amplitude  Hele-Shaw  convection  with  an  imposed  shear.  Their 
interests  were  also  to  determine  the  effect  of  convective  transport  of  momentum  upon 
the  mean  flow.  Prof.  S.  Blumsack  has  continued  his  analysis  of  approximate  solutions 
of  the  barotropic  vorticity  equation  subject  to  periodic  boundary  conditions.  In 
another  study.  Prof.  R.  L.  Pfeffer  and  student  Wen  Ding  have  made  considerable  progress 
solving  the  nonlinear  steady  state  barotropic  vorticity  equation  for  parameters  which 
correspond  to  a  series  of  laboratory  experiments  involving  forced  flow  of  a  rotating 
fludid  over  bottom  topography.  In  these  experiments,  conducted  in  our  Institute,  the 
fluid  is  contained  in  a  rotating  circular  cylindrical  annulus  with  a  differentially 
rotating,  rigid,  radially  sloping,  lid  in  contact  with  the  top  surface  of  the  fluid. 

The  fluid  is  forced  into  motion  by  the  rotation  of  the  lid.  In  another  effort. 
Professors  Barcilon  and  Blumsack  and  Mr.  Y.  Chang  are  investigating  the  dynamics  of 
equilibrium  states  in  a  sheared  barotropic  channel  flow.  This  flow  is  relevant  to 
meteorological  flows  and  the  object  of  this  work  is  to  understand  the  population  of 
certain  states,  called  free  modes  and  comprehend  how  the  actual  flow  in  phase  space 
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new  parallel  computing  techniques  and  fast  algorithms  for  the  efficient  utilization 
of  parallel  computers  in  the  solution  of  the  shallow  water  equations.  Progress  has 
been  made  toward  understanding  the  development  and  nonlinear  behavior  of  traveling 
baroclinic  waves  in  the  atmpsphere  and  their  interaction  with  topographically  and 
thermally  forced  planetary  waves. 
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INTRODUCTION 


The  second  year  of  this  University  Research  Initiative  has  been  even  more  productive  than  the 
first.  To  date,  6  papers  have  been  published  or  accepted  for  publication  and  2  have  been  submitted 
for  publication  based  on  full  or  partial  support  of  AFOSR  89-0462  ais  follows: 


A  study  of  Baroclinic  Wave  Behavior  over  Bottom  Topography  Using  Complex  Principal  Compo¬ 
nent  Analysis  of  Experimental  Data,  by  R.  L.  Pfeffer,  J.  Ahlquist,  R.  J.  Kung,  Y.  Chang  and 
G.  Li;  J.  Atmos.  Sci.,  47,  67-81,  1990. 

•  Effects  of  Wave- Wave  and  Wave- Mean  Flow  Interactions  on  the  Evolution  of  a  Baroclinic  Wave, 

by  A.  Barcilon  and  T.  Nathan;  Geophys.  Astrophys.  Fluid  Dyn.,  56,  59-79,  1991. 

•  Asymmetric  Ekman  Dissipation,  Sloping  Boundaries  and  Linear  Baroclinic  Instability,  by  H. 

Weng  and  A.  Barcilon;  Geophys.  Astrophys.  Fluid  Dyn.,  in  press,  1991. 

•  Reflection  of  Hydrostatic  Mountain  Waves  from  Spatially  Nonuniform  layers,  by  W.  Blumen  and 

A.  Barcilon;  Geophys.  Astrophys.  Fluid  Dyn.,  in  press,  1991. 

•  A  Comparison  of  the  Impact  of  Two  Time-Differencing  Schemes  on  the  NASA/GLAS  Climate 

Model,  by  R.  L.  Pfeffer,  I.  M.  Navon  and  X.  Zou;  .Accepted  for  publication  in  Mon.  Wea. 
Rev.,  1991. 

•  EXSHALL:  A  Turkel-Zwas  Explicit  Large  Time-Step  FORTRAN  Program  for  Solving  the  Shallow- 

Water  Equations  in  Spherical  Coordinates  by  I.  M.  Navon  and  Jian  Yu,  Accepted  for  publica¬ 
tion  in  Computers  in  Geoscience. 

Hele  Shaw  Convection  with  Imposed  Shear,  Part  I,  by  R.  Krishnamurti  and  H.  Yang;  submitted  to 
J.  Fluid  Mech.,  1991. 

Hele  Shaw  Convection  with  Imposed  Shear,  Part  II,  by  H.  Yang  and  R.  Krishnamurti;  to  be 
submitted  to  J.  Fluid  Mech.,  1991. 


The  following  sections  present  a  summary  of  progress  on  research  under  the  grant,  including  the 
work  of  four  Ph.D.  candidates  (one  of  whom  received  his  Ph.D.  degree  within  the  last  year). 

RAYLEIGH-BENARD  CONVECTION 

Thermal  convection  in  a  horizontal  layer  heated  from  below  (“Benard  convection”),  in  the 
turbulent  range,  sometimes  is  accompanied  by  a  large  scale  mean  shear  flow,  so  that  not  only 
heat  but  also  momentum  is  transferred  across  the  layer  (R.  Krishnamurti  and  L.  N.  Howard, 
Proc.  Nat.  Acad.  Sci.  78,1981-1985,  1981).  This  means  shear  flow  appears  to  arise  spontaneously, 
in  a  manner  reminiscent  of  a  bifurcation,  though  it  occurs  when  the  fluid  is  already  in  turbulent 
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motion.  In  a  recent  study,  Prof.  L.  Howard  {Studies  in  Applied  Math.,  83,  273-285, 1990)  obtained 
bounds  on  the  magnitudes  of  the  heat  and  momentum  fluxes  (in  ipfus  of  the  Rayleigh  number, 
which  measures  the  strength  of  the  thermal  driving)  by  an  extension  of  methods  formerly  used  to 
estimate  the  heat  flux  in  the  absence  of  large  scale  flow  (L.  Howard,  J.  Fluid  Medh.,  17,  405-432, 
1963). 

These  last  results  have  now  been  extended  further,  to  the  case  in  which  there  is  also  a  horizontal 
translation  of  one  of  the  boundary  plates  with  respect  to  the  other,  i.e.,  boundary  conditions  which 
in  the  absence  of  any  turbulence  would  be  consistent  with  a  plane  Couette  flow  superimposed 
on  thermal  convection.  Of  course,  a  large  scale  shear  flow  and  associated  momentum  transfer 
is  to  be  expected  when  there  is  mechanical  driving  alone  (turbulent  plane  Couette  flow),  but  its 
occurrence  with  thermal  driving  alone  suggests  that  the  situation  with  both  present  may  also  have 
some  surprises — in  particular  the  momentum  flux  may  be  enhanced  significantly  by  the  presence 
of  convection. 

In  this  case  there  are  two  external  “driving”  parameters,  the  Rayleigh  number  and  a  Reynolds 
number  based  on  the  velocity  of  the  moving  plate.  The  estimate  then  takes  the  form  of  the 
description  of  a  finite  region  in  the  plane  of  heat  flux  and  momentum  flux,  depending  on  the  Rayleigh 
and  Reynolds  numbers,  within  which  any  realized  values  of  these  fluxes  must  lie.  These  results 
have  been  presented  at  the  Miles  Symposium  at  Scripps  Institute  of  Oceanography;  publication  is 
expected  in  the  near  future. 

The  estimates  obtained  depend  on  the  solution  of  a  rather  complicated  variational  problem 
which  has  so  far  not  in  fact  been  exactly  solved,  but  only  itself  bounded  in  terms  of  a  somewhat 
simpler  variational  problem  that  Dr.  Howard  had  solved  before  (loc.  cit.).  It  was  possible  to  solve 
that  problem  analytically  in  terms  of  elliptic  integrals,  but  a  similar  analytical  solution  of  the  more 
complicated  one  has  eluded  him;  a  numerical  study  is  probably  feasible  and  is  under  consideration. 
(Indeed,  a  significant  amount  of  numerical  work  was  required  in  adapting  the  old  elliptic  integral 
solution  to  the  bound  of  the  more  complicated  variational  problem.)  The  difficulty  is  not  so  much 
one  of  a  large  scale  numerical  problem,  but  of  the  need  to  explore  a  larger  parameter  space.  It  is 
perhaps  also  questionable  that  a  marked  improvement  on  the  flux  estimates  will  be  obtained,  even 
if  the  estimate  on  the  variational  problem  is  replaced  by  its  true  solution. 
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In  another  study,  Prof.  R.  Krishnamurti  and  Dr.  H.  Yang  have  completed  their  investigation  of 
finite  amplitude  Hele-Shaw  convection  with  an  imposed  shear.  Their  interests  were  also  to  determine 
the  effect  of  convective  transport  of  momentum  upon  the  mean  flow.  An  imposed  shear  flow  upon 
a  horizontal  layer  of  fluid,  convecting  at  low  Rayleigh  numbers,  leads  to  an  alignment  of  rolls  with 
their  vorticity  vectors  jdong  the  imposed  flow  (“longitudinal  rolls”).  The  Hele-Shaw  geometry,  on 
the  other  hand,  forces  the  important  constraint  that  convective  cells  have  their  vorticity  vectors 
orthogonal  to  the  imposed  Couette  or  Poisieuille  flow  (“transverse  rolls”).  It  is  only  in  this  latter 
arrangement  that  vertical  transport  of  horizontal  momentum  is  found  to  be  counter  to  the  imposed 
gradient.  The  convective  transport  thus  tends  to  enhance  the  imposed  shear. 

The  Hele-Shaw  geometry  has  the  additional  advantage  that  the  flow  is  two-dimensional.  Thus 
Drs.  Yang  and  Krishnamurti  were  able  to  solve  the  Boussinesq  equations  (rather  than  truncated 
model  equations  as  in  Howard  and  Krishnamurti,  1986)  by  using  a  modified  perturbation  expansion 
in  powers  of  the  amplitude  c,  and  a  regular  perturbation  expansion  in  the  powers  of  the  Peclet 
number  representing  the  imposed  shear  flow. 

This  work  is  represented  in  the  Ph.D.  dissertation  of  H.  Yang,  and  two  papers  by  Yang  and 
Krishnamurti,  one  of  which  was  recently  submitted  for  publication. 


LOW  DIMENSIONAL  BAROTROPIC  PROBLEMS 

Prof.  S.  Blumsack  has  continued  his  analysis  of  appro.ximate  solui.ons  of  the  barotropic  vorticity 

equation  -  ^yV-^’x  +  V’l^'V'y  =  0,  in  the  square  0  <  i  <  2k,  0  <  y  <  2k,  subject 

to  periodic  boundary  conditions.  Highly-truncated  doubly  periodic  series  of  the  form  Tp{x,y)  = 

^afc/(t)e'(*®+'*'^  were  assumed,  resulting  in  a  system  of  ordinary  differential  equations  for  the 
k,l 

(complex)  coefficients  aki{t).  This  approach  has  been  undertaken  to  help  check  and  understand 
numerical  models. 

Special  attention  was  devoted  to  a  particular  truncation — a  rhomboidal  truncation  for  which 
1 1^1  “1^1 1  ^  2  and(k,/)  ^  (0,0).  The  resulting  dynamical  system  contains  twelve  degrees  of  freedom. 
A  six  parameter  family  of  equilibrium  states  was  found.  Two  of  the  parameters  account  for  the 
fact  that  any  equilibrium  state  can  be  translated  arbitrarily  in  the  x-y  plane  to  become  another 
equilibrium  state.  The  other  four  parameters  can  be  taken  to  be  the  magnitudes  of  the  Fourier 

/' 
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modes. 


The  time-dependence  for  this  twelve  mode  truncation  has  been  investigated  in  detail.  The 
simplified  dynamics  require  that  the  amplitudes  for  modes  {k,l)  =  (0,  ±2),  (±2,0)  be  constant  in 
time.  Conservation  of  energy  and  enstrophy  then  imply  that  the  total  energy  for  the  subset  of 
modes  {{k,l)  =  (l,±l),(-l,±l)}  is  constant.  In  fact,  this  four-dimensional  subspace  splits  into 
a  pair  of  two-dimensional  subspaces,  with  energy  conserved  in  each.  A  similar  argument  results  in 
energy  conservation  for  the  subset  {(A;,/)  =  (±1,0),(0,  ±1)}.  Trajectories  in  this  four-dimensional 
subspace  are  generally  complicated,  but  several  special  cases  do  result  in  circular  trajectories  in 
certain  planes.  These  cases  have  been  analyzed  in  detail. 

Most  of  the  issues  related  to  the  invariances  and  time  dependence  for  the  twelve  mode  model 
have  been  dealt  with.  The  next  level  model  in  terms  of  comple.xity  has  twenty-four  degrees  of 
freedom;  a  rectangular  truncation  containing  all  modes  such  that  |fc|  <  2,  |/|  <  2  (but  again 
(fc,/)  (0,0)).  Once  again  there  exists  a  six  parameter  family  of  equilibrium  states.  Two  of  these 

are  associated  with  the  translational  invariance  of  the  problem  and  one  of  them  is  associated  with 
the  fact  that  if  V’o  is  an  equilibrium  solution,  then  so  is  any  multiple  of 

Effort  continued  in  the  analysis  of  the  time-dependent  problem.  The  model  was  linearized 
about  an  arbitrary  equilibrium  state  and  the  resulting  matrix  computed.  The  fact  that  the  equilib¬ 
rium  state  can  be  represented  by  means  of  six  parameters  implies  that  six  (out  of  the  twenty-four) 
eigenvalues  must  be  zero.  Eigenvalues  having  non-zero  values  for  both  real  and  imaginary  parts 
come  in  clusters  of  four  of  the  form  ±a±:6;  there  can  be  at  most  four  of  these  clusters  since  there  are 
only  eighteen  non-zero  eigenvalues.  Limited  numericaJ  experimentation  indicated  two  such  clusters 
plus  five  pair  of  complex  conjugate,  purely  imaginary  eigenvalues. 

The  eigenvalue  structure  suggests  that  the  dynamical  system  is  unlikely  to  have  heteroclinic 
orbits  (that  do  not  consist  entirely  of  equilibrium  points).  .\t  most  nine  trajectories  leave  and  at 
most  nine  approach  any  equilibrium  point.  There  seems  to  be  too  much  room  in  our  twenty-four 
dimensional  space  to  permit  a  departing  trajectory  to  meet  with  a  trajectory  incoming  at  a  different 
equilibrium  point. 

In  another  study,  Prof.  R.  L.  Pfeffer  and  student  Wen  Ding  have  made  considerable  progress 
solving  the  nonlinear  steady  state  barotropic  vorticity  equation  for  parameters  which  correspond 


4 


to  a  series  of  laboratory  experiments  involving  forced  flow  of  a  rotating  fluid  over  bottom  topogra¬ 
phy.  In  these  experiments,  conducted  in  our  Institute,  the  fluid  is  contained  in  a  rotating  circular 
cylindrical  annulus  with  a  differentially  rotating,  rigid,  radially  sloping,  lid  in  contact  with  the  top 
surface  of  the  fluid.  The  fluid  is  forced  into  motion  by  the  rotation  of  the  lid.  The  governing 
equation  for  the  steady  state  flow  that  develops  when  we  are  far  from  resonance  is 

J(^,  +  h  -  pR)  +  EV*i>  =  (1) 

T 

top 

surface 
slope 
(/3-effect) 

I 

topographic  Ekman 

forcing  boundary 

layer 

dissipation 

The  side  wall  boundary  conditions  are  d^/dR  =  —  0.  Here,  ^  is  the  dimensionless  stream 

function,  e  is  the  Rossby  number  (ratio  of  inertial  to  Coriolis  acceleration),  h  is  the  ratio  of  the 
topographic  height  to  the  mean  fluid  depth,  /3  is  a  dimensionless  parameter  proportional  to  the 
radial  slope  of  the  rigid  lid  (which  simulates  the  variation  of  the  Coriolis  parameter  with  the  latitude 
on  the  earth),  R  is  radius  from  the  ajcis  of  rotation  and  E  is  the  Ekman  number  (ratio  of  viscous 
force  to  Coriolis  acceleration). 

Dr.  Pfeffer  and  Mr.  Ding  have  compared  solutions  using  these  different  methods  to  determine 
an  acceptable  level  of  truncation.  The  idea  is  to  determine  the  lowest  order  nonlinear  system  that 
can  account  for  the  observed  fluid  behavior.  By  studying  the  solutions  of  such  a  system  they  hope 
to  gain  insight  into  the  physical  mechanisms  involved  in  different  parameter  ranges  and  to  set  the 
stage  for  deriving  a  low  order  dynamical  system  to  account  for  time-dependent  behavior  of  related 
baroclinic  flows.  The  three  methods  used  to  solve  this  problem  are  finite  difference,  semi-spectral 
and  Chebyshev  collocation.  In  the  finite  difference  method  they  solved  the  problem  iteratively  using 
a  grid  of  31  points  in  the  radial  direction  and  61  in  the  azimuthal  direction.  In  the  semi-spectral 
approach  they  used  31  grid  points  in  the  radial  direction  and  5  azimuthal  wave  numbers.  And 
in  the  Chebyshev  collocation  method  they  used  7  Chebyshev  polynomials  to  represent  the  radial 
variation  of  the  stream  function  and  5  Fourier  components  to  represent  the  azimuthal  variation. 


T  T 

internal  forcing  by 
diffusion  differentially 
rotating  lid 


r 

vorticity 

advection 
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Figs,  la,  b  and  c  show  the  solutions  for  the  phase  and  amplitude  of  azimuthal  wave  number  two 
(the  topographic  wave  number)  and  the  azimuthally-averaged  azimuthal  velocity  U,  respectively, 
as  a  function  of  the  rotation  rate  fl  of  the  system.  The  solutions  are  almost  indistinguishable  from 
one  another.  Figs.  Id,  e  and  f  show  experimentally  determined  values  of  the  same  parameters  at 
five  different  rotation  rates  against  a  backdrop  of  the  finite  difference  solution.  The  quantitative 
behavior  of  the  solution  is  similar  to  that  of  the  experimental  data.  There  are,  however,  significant 
quantitative  discrepancies  that  suggest  reexamination  of  certain  assumptions  inherent  in  the  theory 
^  expressed  by  Eq.  (1).  One  such  assumption  is  that  the  effect  of  topography  and  that  of  bottom 
and  top  boundary  layer  dissipation  are  completely  separable  as  they  are  for  small  slopes.  For  large 
slopes,  certain  correction  terms  come  into  play  which  will  have  to  be  taken  into  account.  Another 
assumption  is  that  the  side  wall  boundary  layers  can  be  neglected.  This  becomes  increasingly 
invalid  as  the  rotation  rate  becomes  large.  In  the  coming  year.  Dr.  Pfeffer  and  Mr.  Ding  plan  to 
incorporate  the  appropriate  corrections  in  an  effort  to  bring  the  theory  into  full  compliance  with 
the  experimental  data. 

In  spite  of  the  obvious  quantitative  differences  between  theory  and  experiment,  the  qualitative 
similarities,  and  the  fact  that  they  can  be  represented  well  by  a  low  order  model  in  terms  of  a 
truncated  series  of  Fourier  components  in  the  azimuthal  direction  and  Chebyshev  polynomials  in 
the  radial  direction,  suggests  that  it  should  be  possible  to  use  the  low  order  model  to  shed  light  on 
the  relative  roles  of  the  /J-effect  and  dissipation  in  bringing  about  the  variations  of  phase,  amplitude 
and  U  with  Q  shown  in  Fig.  1.  This  will  be  the  focus  of  Dr.  Pfeifer’s  and  Mr.  Ding’s  efforts  in  the 
coming  year. 

Some  questions  that  have  arisen  in  connection  with  Dr.  Pfeifer’s  and  Mr.  Ding’s  work  on  this 
problem  have  led  Prof.  Howard  to  look  more  closely  into  the  derivation  of  approximate  equations 
and  their  boundary  conditions  for  rotating,  nearly  geostrophic,  flow  in  containers  of  arbitrary  shape 
and  with  rather  general  boundary  conditions.  He  had  previously  considered  such  questions  in  some 
detail  (also  for  stratified  flow  (L.  N.  Howard,  Lectures  in  Applied  Math.,  13,  121-124,  1971),  but 
only  in  the  linear  context.  His  conjecture  is  that  the  non-linear  quasi-geostrophic  system,  though 
it  has  been  extensively  studied  and  applied,  has  not  yet  been  quite  properly  fitted  into  a  general 
formulation.  This  investigation  is  not  finished  yet,  but  has  reached  a  stage  where  -  preliminary 
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Fig.  1.  Solutions  of  the  barotropic  vorticity  equation  by  three  different  methods  for  a.  phase 
vs  n,  b.  amplitude  vs  0.  and  c.  U  vs  il.  The  slight  discrepancies  disappear  completely  when  we 
take  more  grid  points  in  the  azimuthal  direction  in  the  finite  difference  solution. 

Comparisons  between  the  finite  difference  solution  and  e.xperimental  determinations  of  d. 
phase,  e.  amplitude  and  c.  azimuthal  mean  velocity  U. 


attempt  at  writing  it  up  seems  desirable,  and  Dr.  Howard  plans  to  do  this  in  the  coming  year.  The 
basic  structure  of  the  formulation  can  be  described  briefly  as  follows. 

When  the  Rossby  and  Ekman  numbers  are  small,  and  time  scales  are  long  compared  to  the 
rotation  period,  the  interioi  flow  approximately  satisfies  the  (homogeneous)  linear  geostrophic  equa¬ 
tions,  with  approximately  zero  normal  velocity.  However,  the  homogeneous  geostrophic  equations, 
with  zero  normal  velocity,  have  in  general  many  solutions;  in  the  simplest  case  where  the  “depth” 
of  the  container — measured  parallel  to  the  rotation  axis — is  constant,  the  pressure  is  essentially 
an  arbitrary  function  of  the  transverse  variables  (and  determines  the  velocity  components).  When 
the  depth  is  not  constant,  the  arbitrariness  is  less  extreme  (the  geostrophic  flow  must  be  along 
contours  of  constant  depth),  but  still  present.  As  will  all  linear  problems,  when  the  homogeneous 
form  has  many  solutions,  the  non-homogeneous  equations  (and  boundary  conditions)  can  have  so¬ 
lutions  only  when  the  non-  homogeneities  satisfy  certain  “solvability  conditions."  These  conditions 
are  sometimes  expressed  as  conditions  of  orthogonality  to  all  solutions  of  the  adjoint  homogeneous 
problem;  in  the  present  case  the  adjoint  homogeneous  problem  has  many  solutions,  depending  in 
part  on  the  nature  of  the  container,  and  the  solvability  conditions  are  somewhat  complicated  in 
general.  The  full  equations  can  be  regarded  as  the  geostrophic  equations  with  non-homogeneities 
given  by  the  terms  involving  the  Rossby  and  Ekman  numbers.  The  boundary  conditions  on  the 
interior  flow  come  from  the  imposed  boundary  conditions  (differential  motions  of  the  boundary, 
sources  and  sinks,  etc.)  modified  by  the  effects  of  boundary  layers  (whose  structure  is  in  general 
also  coupled  to  the  interior  flow)  and  also  appear  as  inhomogeneities  in  the  geostrophic  problem. 
The  solvability  conditions  then  lead  to  additional  relations  on  the  interior  flow  (like  the  potential 
vorticity  equation  and  the  Ekman  layer  matching — these  relations  are  not  necessarily  linear),  which 
ultimately  determine  which  of  the  many  geostrophic  flows  is  in  fact  the  first  approximation  to  the 
reaJ  flow,  and  what  its  higher  corrections  are.  This  general  description  is  fairly  straightforward, 
but  working  it  out  in  detail  is  rather  complex. 

In  another  effort.  Professors  Barcilon  and  Blumsack  and  Mr.  Y.  Chang  are  investigating  the 
dynamics  of  equilibrium  states  in  a  sheared  barotropic  channel  flow.  This  flow  is  relevant  to 
meteorological  flows  and  the  object  of  this  work  is  to  understand  the  population  of  certain  states, 
called  free  modes  (see  below),  and  comprehend  how  the  actual  flow  in  phase  space  is  attracted  from 
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one  state  to  another.  Thus,  this  problem  has  some  bearing  on  atmospheric  predictability  and  was 
strongly  motivated  by  the  work  of  Branstator  and  Opsteegh  ( J.  Atmos.  Set.,  46,  1799-1814,  1989). 

Thus,  they  are  interested  in  finding  the  free  states  of  the  barotropic  vorticity  equation  for  a 
given  shear  flow.  These  states  are  characterized  by  the  equation 

=  0 

where  ^  is  the  streamfunction  and  V'*}  the  vorticity. 

Following  VaJlis  et  al.  {J.  Fluid  Slech.,  207,  133-152,  1989)  they  have  developed  a  numerical 
scheme  which  locates  these  states  starting  from  a  given  initial  condition.  The  property  of  this 
scheme  is  that  it  conserves  the  vorticity  of  a  given  particle  but  allows  the  total  energy  of  the  system 
to  be  a  monotonic  function  of  time,  the  total  energy  approaching  a  constant  when  a  free  state  is 
reached. 

Research  has  been  pursued  along  two  avenues: 

1.  They  determined  possible  steady  wavy  solutions  of  the  barotropic  vorticity  equation.  These 
solutions  involve  sheared  flows  constructed  by  means  of  trigonometric  functions. 

2.  They  chose  as  an  initial  basic  state 


“0(2/)  =  -«o  tanh(y) 

superimposed  an  initial  perturbation  of  the  form 

i'{x,y,Q)  -  .Aosinioiy  +  n)s'm koX 

and  integrated  the  vorticity  equation  using  a  spectral  method.  Various  initial  conditions  were 
considered  by  specifying  different  values  of  the  zonal  wave  number  ko  (the  most  unstable  ko 
=  0.44  by  linear  theory)  as  well  cis  different  numbers  of  modes.  They  then  compared  these 
dynamics  with  the  dynamics  of  the  packages  that  search  for  the  free  states  of  the  system. 

A  short  note  entitled  “Some  Preliminary  Results  on  Equilibrium  States  in  a  Sheared  Barotropic 
Flow”  by  Y.  Chang,  A.  Barcilon  and  S.  Blumsack  was  sent  earlier  to  AFOSR.  The  authors  plan  to 
expand  this  note  into  a  paper  which  will  be  submitted  for  publication  before  the  end  of  the  year. 
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NUMERICAL  METHODS 


The  shallow  water  equations  are  a  set  of  first  order  nonlinear  hyperbolic  parti  J  differential 
equations  having  many  important  applications  in  meteorology  and  oceanography.  These  equations 
can  be  used  in  studies  of  tides  and  surface  water  run-off.  They  may  also  be  used  to  study  large- 
scale  waves  in  the  atmosphere  and  ocean  if  terms  representing  the  effects  of  the  earth’s  rotation 
(Coriolis  terms)  are  included.  Indeed  it  has  become  customary,  in  developing  new  numerical  meth¬ 
ods  for  numerical  weather  prediction  or  oceanography,  to  study  first  the  simpler  nonlinear  shallow 
water  equations  system,  which  possesses  the  same  mixture  of  slow  and  fast-moving  waves  as  the 
more  complex  baroclinic  three-dimensional  primitive  equations  of  motion.  Galerkin  finite  element 
methods  have  been  successfully  applied  to  these  equations  by  many  researchers.  However,  these 
algorithms  were  not  designed  to  run  efficiently  on  parallel  processing  architectures. 

In  recent  years  there  was  a  great  interest  in  parallel  computations  due  to  the  advent  and  growing 
popularity  of  parallel  computers.  With  an  appropriate  design  of  the  algorithm,  all  the  processors 
within  a  single  computer  can  be  called  upon  to  work  on  a  single  program  concurrently  and  thus 
the  result  can  be  obtained  in  a  much  shorter  period  of  time.  A  parallel  computer  with  n  processors 
is  cheaper  in  CPU  time  than  n  computers  each  having  a  single  processor.  Parallel  computers  can 
be  used  to  tackle  such  problems  ais  worldwide  weather  forecasting  and  the  like  which  constitute  a 
complex  ta^k  for  serial  machines.  By  carrying  out  research  in  parallel  processing,  research  workers 
may  further  explore  the  potential  power  of  computer  resources.  Parallel  processing  makes  it  possible 
to  attain  the  goal  of  full  utilization  of  the  most  powerful  high-performance  supercomputers. 

In  view  of  the  importance  of  the  shallow  water  equations  in  meteo.ology  and  oceanography  and 
the  obvious  advantages  offered  by  modern  parallel  computers,  the  goal  in  the  present  research  effort 
is  to-investigate  new  parallel  computing  techniques  and  develop  fast  algorithms  for  the  efficient 
utilization  of  parallel  computers  in  the  solution  of  the  aforementioned  set  of  partial  differential 
equations.  The  expected  results  of  this  research  will  impact  on  both  the  area  of  parallel  scientific 
computation  as  well  as  on  application  problems  in  meteorology  and  oceanography. 

Parallelism  present  in  the  numerical  solution  algorithms  of  PDE’s  can  be  exploited  at  several 
levels.  Domain  decomposition  is  the  highest  level  of  parallelism  in  the  sense  that  the  whole  problem 
is  split  up  into  several  smaller  sub-problems  which  will  be  apportioned  among,  and  handled  by. 
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several  processors  concurrently  even  before  we  discretize  it  to  obtain  the  appropriate  linear  or  non¬ 
linear  systems  of  algebraic  equations.  This  technique  has  been  chosen  as  the  general  methodology 
of  Prof.  I.  M.  Navon  and  student  Yihong  Cai  to  be  employed  in  mapping  the  whole  computational 
effort  into  a  number  of  processors  for  parallel  processing. 

Under  the  present  URI,  they  have  successfully  decomposed  the  limited  area  domain,  over 
which  the  finite-element  model  of  the  hyperbolic  nonlinear  shal’ow  water  equations  are  defined, 
into  two  subdomains  and  then  into  four  subdomains  by  applying  a  Schur  complement  domain 
decomposition  technique,  which  is  a  nonoverlapping  domain  decomposition  approach.  The  results 
were  more  satisfactory  than  expected  in  the  sense  that  a  speed-up  had  already  been  achieved  even 
on  a  single  processor.  They  have  also  formulated  the  framework  of  the  general  technique  which  can 
be  applied  to  an  arbitrary  number  of  subdomains. 

One  of  the  key  steps  for  carrying  out  the  Schur  domain  decomposition  is  to  obtain  the  arrow¬ 
head  shaped  matrix,  from  which  information  can  be  derived  regarding  the  Schur  matrix  complement 
and  the  systems  governing  the  various  subdomain  problems.  To  build  a  device  into  their  program 
which  allows  for  the  easy  modification  of  the  code  to  deal  with  different  resolutions  for  obtaining 
higher  ana  higher  orders  of  accuracy  and  testing  for  the  convergence  of  the  method.  Dr.  Navon  and 
Mr.  Cai  have  found  and  devised  an  efficient  way,  namely  a  renumbering  algorithm,  for  transforming 
the  original  system  matrix  into  the  arrow-head  matrix.  This  renumbering  algorithm  allows  for  an 
easier  implementation  designed  to  carry  out  the  Schur  complement  domain  decomposition  and  also 
provides  an  easy  means  by  which  to  efticiently  implement  multicoloring  techniques. 

Another  very  important  issue  in  the  Schur  complement  domain  decomposition  is  how  to  avoid 
the  explicit  evaluation  of  the  Schur  complement.  It  turns  out  that  is  very  e.xpensive  to  obtain  the 
Schur  complement  as  it  requires  kn,  solves,  where  n,  is  the  number  of  nodes  in  the  interface  and 
k  is  the  number  of  subdomains. 

The  explicit  computation  of  the  Schur  complement  will  probably  dominate  the  overall  solution 
process  for  large  degrees  of  freedom  on  the  interface  as  the  mesh  is  refined.  To  apply  such  accelera¬ 
tion  methods  as  Conjugate  Gradient  Squared  (CGS)  or  Generalized  Minimal  Residual  (GMRES)  to 
the  nonsymmetric  matrix  interface  system.  Dr.  Navon  and  Mr.  Cai  have  formulated  an  algorithm 
which  avoids  the  e.xplicit  calculation  of  Schur  complement.  This  algorithm  can  be  described  as  the 
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inner  loop  across  the  subdomains  and  outer  loop  on  the  interface.  They  are  now  at  the  stage  where 
they  are  ready  to  carry  out  extensive  numerical  experiments  based  on  this  idea  and  the  results  will 
be  published  in  the  coming  months. 

The  interface  system  is  usually  solved  by  an  iterative  technique  and  each  iteration  for  the 
interface  requires  solving  all  the  subdomain  problems  once.  It  is  therefore  of  paramount  importance 
to  make  the  number  of  iterations  on  the  interface  as  small  as  possible.  In  order  to  achieve  this, 
various  interface  preconditioning  techniques  will  be  investigated  in  the  coming  year  so  as  to  change 
the  eigen-spectrum  properties  in  order  to  reduce  the  condition  number  of  the  problem  as  well  as  to 
cluster  the  eigenvalues  in  the  spectrum  and  thus  speed  up  the  iteration  process. 

Both  modified  incomplete  LU  decomposition  preconditioning  methods  as  well  as  boundary 
probe  methods  will  be  considered  as  prime  candidates.  This  work  is  important  since,  if  an  optimal 
preconditioner  can  be  found  for  either  the  interface  system  or  the  subdomain  nonsymmetric  matrix 
systems,  the  time  required  for  calculation  will  be  dramatically  reduced.  This  is  because  there  will 
result  a  significant  reduction  in  the  number  of  iterations  required  to  achieve  a  pre-specified  conver¬ 
gence  criterion.  Since  the  solution  of  the  interface  system  is  the  overhead  for  parallel  processing, 
they  will  also  study  the  possibilities  of  parallelizing  both  the  acceleration  methods  as  well  as  the 
preconditioners.  A  further  speed-up  can  be  expected  from  this  research. 

An  additional  effort  will  be  dedicated  towards  parallelizing  the  various  phases  in  the  finite 
element  discretization  algorithm  such  as  the  assembly  procedure,  node  numbering  and  calculation 
of  the  element  matrices.  The  computational  advantage  of  exploiting  this  type  of  parallelism  has 
already  attracted  the  attention  of  many  researchers.  However,  there  is  still  a  sizable  amount  of 
research  to  be  carried  out  for  obtaining  adequate  parallel  algorithms.  These  parallel  algorithms 
will  be  tested  on  various  computer  architectures  such  as  the  Cray  Y-MP  supercomputer,  hypercube 
machines  and  Alliant  parallel  computers. 

So  far,  Dr.  Na^'on  and  Mr.  Cai  we  have  carried  out  the  domain  decomposition  using  the  nonover¬ 
lapping  approach.  However,  there  is  another  domain  decomposition  approach  called  the  overlapping 
approach.  This  consists  in  a  broad  class  of  techniques  in  which  initial  guesses  are  required  for  the 
overlapped  subdomain  regions.  The  solution  involves  an  iterative  process  to  be  carried  out  until 
the  difference  between  different  subdomains  over  the  same  overlapped  region  achieves  a  prescribed 
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tolerance.  They  plan  to  use  the  theory  of  optimal  control  to  address  this  problem  in  the  future. 
This  will  require  them  to  develop  the  adjoint  model  of  the  finite-element  discrete  version  of  the 
shallow-water  equations  and  to  use  a  large  scale  unconstrained  minimization  algorithm  to  mini¬ 
mize  a  least  squares  cost  functional  measuring  the  lack  of  fit  between  two  overlapped  subdomain 
solutions. 


BAROCLINIC  INSTABILITY 

This  section  deals  with  our  progress  toward  understanding  the  development  and  nonlinear 
behavior  of  traveling  baroclinic  waves  in  the  atmosphere  and  their  interaction  with  topographically 
and  thermally  forced  planetary  waves. 

In  a  paper  entitled  “Effects  of  Wave- Wave  and  Wave-Mean  Flow  Interactions  on  the  Evolution 
of  a  Baroclinic  Wave”  (Geophys.  Astrophys.  Fluid  Dyn.,  56,  59-79,  1991),  Professors  A.  Barcilon 
and  T.  Nathan  (visitor  from  the  University  of  California  at  Davis)  examined  within  the  context 
of  a  two-layer  model  the  weakly  nonlinear  evolution  of  a  free  baroclinic  wave  in  the  presence  of 
slightly  supercritical,  vertically  sheared  zonal  flow  and  a  forced  stationary  wave  field  that  consists 
of  a  single  zonal  scale  and  an  arbitrary  number  of  meridional  harmonics.  The  presence  of  the 
(planetary-sciile)  stationary  wave  introduces  zonal  variations  in  the  supercriticality  and  is  shown 
to  alter  the  growth  rate  and  asymptotic  equilibrium  of  the  (synoptic-scale)  baroclinic  wave  via 
two  distinct  mechanisms:  The  first  is  due  to  the  direct  interaction  of  the  stationary  wave  with 
that  portion  of  the  mean  field  that  is  corrected  by  the  zonally  rectified  stationary  wave  fluxes 
(wave-mean  mechanism).  These  mechanisms  can  oppose  or  augment  each  other  depending  on  the 
amplitude  and  spatial  structure  of  the  stationary  wave  field.  If  the  stationary  wave  field  is  confined 
primarily  to  the  upper  (lower)  layer  and  consists  only  of  the  gravest  cross-stream  mode,  conditions 
are  favorable  (unfavorable),  for  nonzero  equilibrium  of  the  free  wave. 

In  addition  to  the  time-dependent  heat  flux  generated  by  baroclinic  growth  of  the  free  wave,  its 
interaction  with  a  stationary  wave  field  consisting  of  two  or  more  meridional  harmonics  generates 
time-dependent  heat  fluxes  that  vary  with  the  period  of  the  free  wave.  However,  if  the  stationary 
wave  field  contains  several  meridional  harmonics  of  sufficiently  large  amplitude,  the  free  baroclinic 
wave  is  destroyed. 
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In  another  paper  entitled  “Asymmetric  Ekman  Dissipation,  Sloping  Boundaries  and  Linear 
Baroclinic  Instability”  {Geophys.  Astrophys.  Fluid  Dyn.,  54,  1991,  in  press)  Drs.  H.-Y.  Weng  and 
A.  Barcilon  examined  the  effects  of  symmetric  and  asymmetric  Ekman  dissipation  on  baroclinic 
instability,  phase  speed  and  wave  structure  in  a  linear  Eady-like  channel  model  with/ without  oppo¬ 
sitely  sloping  boundaries.  In  the  absence  of  sloping  boundaries,  symmetric  Ekman  dissipation  has 
a  stabilizing  tendency  for  all  waves.  When  the  Ekman  layers  are  asymmetric,  the  viscous  asym¬ 
metry  destabilizes  waves  by  (i)  extending  the  unstable  waveband  toward  both  the  long  and  short 
waves,  and  (ii)  increasing  their  growth  rates,  compared  with  the  viscous  symmetric  Ccise.  When  the 
asymmetric  dissipation  is  very  small,  the  destabilization  may  result  in  a  frictional  instability  for 
short  waves  which  are  stable  in  the  inviscid  case.  The  viscous  asymmetry  makes  waves  dispersive 
and  their  structure  asymmetric  about  mid-depth.  However,  the  sense  of  the  viscous  asymmetry 
and  the  sign  of  shear  do  not  affect  the  instability,  but  do  modify  the  direction  of  phase  propagation 
and  the  shape  of  the  wave  structure. 

With  asymmetric  Ekman  dissipation  and  sloping  boundaries,  frictional  instability  is  a  mixture 
of  three  mechanisms:  (i)  symmetric  dissipation  in  the  presence  of  sloping  boundaries,  (ii)  viscous 
asymmetry  in  the  absence  of  sloping  boundaries;  and  (iii)  the  sense  of  viscous  asymmetry  in  the 
presence  of  sloping  boundaries,  which  is  sensitive  to  the  wavenumber  and  the  sign  of  shear.  The 
slope  renders  the  phase  speed  more  dispersive  and  sensitive  to  the  sign  of  shear  and  the  sense  of 
the  viscous  asymmetry.  For  westerly  shear,  the  slope  reduces  (increases)  the  wave  amplitude  and 
makes  the  wave  more  baroclinic  (barotropic)  in  the  lower  (upper)  level.  For  a  given  wavenumber, 
easterly  shear  dynamics  may  be  deduced  from  westerly  shear  dynamics  by  a  proper  change  of  the 
sense  of  the  viscous  asymmetry. 

In  another  ongoing  project  under  the  URI,  one  of  our  regular  visitors.  Prof.  T.  Nathan  of  U. 
C.  Davis,  is  examining  the  response  of  a  multiwave  baroclinic  wave  system  to  periodic  (seasonal) 
variations  in  zonal  mean  thermal  forcing.  The  goal  is  to  understand  the  role  of  nonlinearity  in  the 
non-cyclic  response  of  the  earth’s  atmosphere  (which  exhibits  strong  differences  in  circulation  and 
climate  from  year  to  year)  to  periodic  forcing  by  the  sun.  In  view  of  the  fact  that  even  the  simplest 
dynamical  system  can  exhibit  a  wide  range  of  spatial  and  temporal  variability  when  subjected 
to  periodic  forcing,  the  following  question  naturally  arises:  to  what  extent  are  the  pronounced 
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differences  in  weather  phenomena  of  the  same  season  of  different  years  due  to  secisonal  forcing? 
Addressing  this  question  forms  the  basis  of  Dr.  Nathan’s  research  under  the  URL 

The  model  is  basically  that  of  Nathan  (J.  Atmos.  Set.,  45,  1052-1071,  1988;  J.  .Atmos.  Sci., 
45,  113-130,  1989),  modified  to  account  for  periodic  variations  in  zonal  mean  potential  vorticity. 
Briefly,  he  considers  the  classic  two-layer  model  for  flow  confined  to  a  zonally  periodic,  midlatitude 
beta-plane  channel.  The  flow  is  forced  by  sinusoidal  bottom  topography  and  periodic  variations  in 
zonal  mean  potential  vorticity.  Potential  vorticity  damping  is  the  sole  dissipative  mechanism. 

The  basic  state  flow  is  zonal,  steady,  and  independent  of  latitudinal  variations.  Dr.  Nathan 
considers  a  disturbance  field  consisting  of  a  free  marginally  unstable  baroclinic  wave  and  a  reso¬ 
nant  topographic  wave.  The  analysis  hinges  on  introducing  “long”  time  scales  and  expanding  the 
disturbance  field  in  an  asymptotic  series.  The  problem  then  reduces  to  a  sequence  of  inhomoge¬ 
neous  linear  differential  equations  where  the  nonlinear  terms  appear  ais  forcing  terms  on  the  linear 
operator.  To  insure  the  validity  of  the  asymptotic  expansions,  one  must  remove  those  terms  which 
are  orthogonal  to  the  homogeneous  adjoint  solution.  This  solvability  condition  yields  the  sought 
after  temporal  evolution  equations  for  the  wave  amplitudes  and  mean  field  corrections. 

The  equations  governing  the  temporal  evolution  of  the  free  unstable  wave  amplitude,  A(f), 
forced  wave  amplitude,  B{t)  and  mean  field  corrections,  $i(y,  f)  and  have  been  derived 

and  can  be  written  in  general  form  as: 
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+  h  ^^6  ^  sin-  iry  +  *  +  H  sin(</r);  j  =  1, 2. 


where, 

M,  =  ^  +  (-1)>F($i-$2)  (4) 

is  the  zonal  mean  potential  vorticity  correction.  In  the  above  equations  Op  (p  =  1-19),  6^  = 

1-9),  and  €„  (n  =  1-6)  are  constraints  that  depend  on  zonal  wavenumbers  of  the  free  and  forced 
waves.  Internal  rotational  Froude  number,  F,  planetary  vorticity  factor,  0,  and  basic  state  zonal 
wind  in  the  upper  layer,  Ui.  The  topographic  height  is  given  by  h,  r  measures  the  strength  of 
the  potential  vorticity  damping,  the  asterisk  represents  a  comple.x  conjugation,  i  =  H  is  the 
amplitude  of  the  zonal  mean  periodic  forcing,  and  r  is  the  period  of  the  zonal  mean  forcing. 

To  cast  Eqs.  (l)-(3)  in  a  form  suitable  for  numerical  integration,  the  mean  field  corrections 
were  written  as 

St 

cosm-y;  ;  =  1,2.  (5) 

m=l 

where  the  terms  with  m  odd  are  zero.  Insertion  of  (5)  into  (1)-  (4)  yields  a  system  of  coupled, 
nonlinear,  ordinary  differential  equations,  which  were  integrated  forward  in  time  using  a  fourth- 
order  scheme  with  a  nondimensional  time  step  of  0.1.  For  the  majority  of  integrations,  Nt  =  12 
proved  sufficient.  The  following  initial  conditions  were  used:  B  =  dAjdt  =  0  at  t  =  0. 

In  the  absence  of  periodic  (seasonal)  forcing,  numerical  integrations  of  the  evolution  equations 
performed  under  the  URI  thus  far  reveal  two  distinct  asymptotic  states  of  the  system:  one  is 
a  single  (stationary)  topographic  wave  state,  or  a  mixed  wave  state  where  the  baroclinic  wave 
propagates  at  a  fixed  amplitude  while  the  topographic  wave  remains  stationary  and  phase  shifted 
with  respect  to  the  topography.  The  single  topographic  wave  states  are,  relatively  speaking,  favored 
for  large  zonal  scales,  large  topographic  heights  and  weak  dissipation.  However,  for  sufficiently  weak 
dissipation,  inclusion  of  zonal  mean  periodic  forcing  produces  quasi-periodic,  mixed  wave  regimes 
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which  correspond  to  “interannual”  variations  in  the  model  circulation.  Wave- wave  interactions  have 
been  shown  to  play  a  vital  role  in  determining  the  “interannual”  variability  of  the  model  circulation. 
In  particular,  wave-wave  interactions  are  responsible  for  the  following:  i)  the  interchange  of  wave 
dominance  from  one  “winter”  (“summer”)  to  the  next  (when  the  system  exhibits  nonperiodic 
responses);  ii)  reducing  some  sub-cycle  responses  to  a  periodic  response  with  the  period  of  the 
forcing;  and  iii)  confining  the  topographic  wave  to  a  more  narrow  region  of  phase  space  than  would 
otherwise  occur  in  the  absence  of  wave-wave  interactions. 

During  the  next  year  Dr.  Nathan  plans  to  extend  the  above  model  to  include  zonally  zisym- 
metric  periodic  forcing  and  simulate  the  influence  of  seasonal  changes  in  land-sea  heating  contrasts. 

MOUNTAIN  WAVES 

Other  research  which  has  developed  in  the  course  of  the  present  URI  has  led  to  an  additional 
paper  entitled  “Reflection  of  Hydrostatic  Mountain  Waves  from  Spatially  Nonuniform  Layers”  ( Geo- 
phys.  Astrophys.  Fluid  Dyn.,  1991,  in  press).  In  this  work.  Professor  A.  Barcilon  has  collaborated 
with  Prof.  W.  Blumen  of  the  University  of  Colorado  in  an  investigation  of  the  steady,  hydrostatic, 
inviscid,  Boussinesq  flow  of  a  stably  stratified  fluid  over  a  bell-shaped  ridge  within  the  framework 
of  a  linear  model.  The  three  layer  model  atmosphere  introduced  is  such  that  the  Brunt-Vaisala 
frequency  is  constant  in  each  layer  but  the  interfaces  of  the  middle  layer  are  allowed  to  vary  gently 
in  the  cross-ridge  direction.  In  essence,  the  model  can  be  tuned  in  both  vertical  and  horizontal 
directions.  These  cross-ridge  variations  can  produce  significant  differences  in  both  the  cross-ridge 
surface  wind  and  the  surface  drag  compared  to  the  response  obtained  by  use  of  a  horizontcilly  uni¬ 
form  reflecting  layer.  These  changes  are  sensitive  to  both  the  vertical  location  of  the  middle  layer 
and  to  the  slope  of  its  lower  interface  at  the  ridge  crest.  Many  of  these  features  are  e.xplained  by 
means  of  a  conventional  layered- model  analysis. 
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